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ABSTRACT 

We use a large set of cosmological smoothed particle hydrodynamics (SPH) sim- 
ulations to examine the effect of mass resolution and box size on synthetic Lya forest 
spectra at 2 < z < 5. The mass resolution requirements for the convergence of the 
mean Lya flux and flux power spectrum at z = 5 are significantly stricter than at 
lower redshift. This is because transmission in the high redshift Lya forest is pri- 
marily due to underdense regions in the intergalactic medium (IGM), and these are 
less well resolved compared to the moderately overdense regions which dominate the 
Lya forest opacity at z ~ 2 — 3. We further find that the gas density distribution 
in our simulations differs significantly from previous results in the literature at large 
overdensities (A > fO). We conclude that studies of the Lya forest at z = 5 using 
SPH simulations require a gas particle mass of M gas < 2 x f0 5 Mo/h, which is >8 

times the value required at z — 2. A box size of at least 40Mpc/h is preferable at all 
redshifts. 

Key words: methods: numerical - intergalactic medium - quasars: absorption lines. 



1 INTRODUCTION 

Most observations of the Lya forest are at 2 < z < 4; 
this is where the highest quality optical spectra and the 
largest data sets are available. Consequently, many stud- 
ies of hydrodynamical Lya forest simu lations also fo cus on 
this redshift range (for a review see Mciksin 2007). The 
convergence of a variety of simulated Lya forest statis- 
tics with resolution and box size has been explored in 
detail at z < 4, both for Eule r ian g r id and Lagrangian 
SPH simulations l|Theuns et all Il998l ; iBrvan et all 1 19991 : 
iRegan et al.ll2007l ). In recent years, however, the importance 
of the z > 4 L ya forest as a p r obe of the hydroge n reioni- 
sation epoch i|Fan et alj l200d : iBecker et al1l2007l ) has led 
to considerable interest in the properties of hydrodynamical 
Lya forest simulations at the highest observable r edshifts 
i|Paschos fc Norman|[2005l : iBolton fc Haehnelu[2007l ). 

However, it is not obvious that the numerical require- 
ments for simulating the Lya forest at z > 4 are the same as 
at lower redshifts. The average transmission of the Lya for- 
est decreases towards higher redshift as the physical gas den- 
sity increases and the intensity of the ultraviolet (UV) back- 
ground falls. Absorption lines associated with the mildly 
overdense regions of the IGM at z ~ 2 are gradually replaced 
by transmission peaks from underdense regions, culminating 
in completely saturated absorption by z ~ 6. The decrease in 



the characteristic overdensity associated with the Lya forest 
towards higher redshift will place strong demands on SPH 
simulations, which naturally focus on resolving high density 
regions, as a function of time. 

In this letter we analyse a large set of cosmological simu- 
lations perfo rmed with an u pgraded version of the SPH code 
GADGET-2 (|Springelll2005l ). We consider the convergence of 
two widely used Lya forest flux statistics, the mean flux 
and the flux power spectrum, with mass resolution and box 
size. We also examine the gas density distribution in our 
simulations in detail. Although some of the issues discussed 
in this letter are generally appreciated, there have been no 
studies of SPH Lya forest simulations at z ~ 5 where they 
are highlighted explicitly. Our intention is that these results 
will provide a useful reference for future modelling of the 
Lya forest at the highest observable redshifts. 

2 SIMULATIONS 

We perform 24 hydrodynamical simulations using GADGET- 
3, an upgraded version of the publicly avai lable paral- 
lel Tree-SPH code GADGET-2 i|Springell [20051 ). The sim- 
ulations are summarised in Table [TJ and cover a wide 
range of comoving box sizes and gas particle masses. 
The simulations are all started at z = 99, with ini- 
tial conditions constructed using the same random seed. 
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Table 1. Mass resolution and box size (comoving) of the hydro- 
dynamical simulations used in this work. 



Name 


L 


Total particle 


iWgas 




[Mpc/h] 


number 


[M /h] 


2.5-50 


2.5 


2 X 50 3 


1.61 x 10 6 


2.5-100 


2.5 


2 X 100 3 


2.01 x 10 s 


2.5-200 


2.5 


2 X 200 3 


2.51 x 10 4 


2.5-400 


2.5 


2 X 400 3 


3.14 x 10 3 


5-50 


5 


2 x 50 3 


1.29 x 10 7 


5-100 


5 


2 x 100 3 


1.61 x 10 6 


5-200 


5 


2 x 200 3 


2.01 x 10 s 


5-400 


5 


2 x 400 3 


2.51 x 10 4 


10-50 


10 


2 x 50 3 


1.03 x 10 s 


10-100 


10 


2 x 100 3 


1.29 x 10 7 


10-200 


10 


2 x 200 3 


1.61 x 10 6 


10-400 


10 


2 x 400 3 


2.01 x 10 5 


20-50 


20 


2 x 50 3 


8.22 x 10 s 


20-100 


20 


2 x 100 3 


1.03 x 10 s 


20-200 


20 


2 x 200 3 


1 29 x 10 7 


20-400 


20 


2 x 400 3 


1.61 x 10 6 


40-50 


10 


2 x 50 3 


6.58 x 10 9 


40-100 


40 


2 x 100 3 


8.22 x 10 s 


40-200 


40 


2 x 200 3 


1.03 x 10 s 


40-400 


40 


2 x 400 3 


1.29 x 10 7 


80-50 


SO 


2 x 50 3 


5.26 x 10 10 


80-100 


80 


2 x 100 3 


6.58 x 10 9 


80-200 


80 


2 x 200 3 


8.22 x 10 s 


80-400 


80 


2 x 400 3 


1.03 x 10 s 



The cosmological parameters are (f2 m , Qa, ftbh 2 , h, erg, n B ) = 
(0.26, 0.74, 0.024, 0.72, 0.85, 0.95). These are consistent with 
the fifth year Wilkinson Microw ave Anisotropy Probe 
(WMAP) data l|Dunklev et al.ll2009T l aside from a s, which is 
in better agreement wi th Lya forest constraints (|Viel et al.l 
l2004l ; ISeliak et al.ll2005l ). The gas is assumed to be of primor- 
dial composition with a helium mass fraction of Y = 0.24. 
The gravitational softening length is set to l/30 th of the 
mean linear interparticle spacing. Star formation is included 
using a simplified prescription which converts all gas parti- 
cles with overdensity A = p/(p) > 10 3 and temperature 
T < 10 5 K into collisionless stars. 

The baryons in the simulations are photoionised and 
heated by the UV background model of lHaardt fc Mad"aul 
|200ll ) which includes emission from both quasars and galax- 
ies. The UV background is switched on at z = 9 and is ap- 
plied in the optically thin limit using a non-equilibrium ion- 
isation algorithm. The He n photo- heating rate is increased 
by a factor of 1.8 to give temperatu res similar to exist- 
ing measurements (jSchave et alj 1200(3) . Synthetic Lya for- 
est spectra are constructed from each simulation at z = 
(2,3,4,5). Note that in this work we do not subsequently 
alter the resolution or S/N of the spectra. 

3 RESULTS 

3.1 Convergence of flux statistics 

We first consider the simplest Lya forest observable, the 
mean flux, (F) — (e _T ), in Fig. [T] This shows the differ- 
ence in (F) measured from our synthetic spectra relative to 
the highest resolution model (2.5-400) as a function of gas 
particle mass at z = (2,3,4,5). At z = 2 the mean flux 
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Figure 1. The difference in the mean Lya flux, SF = ((F) — 
(F){ )/(F)f relative to a fiducial value, (F){, as a function of gas 
particle mass for the hydrodynamical simulations listed in Ta- 
ble [l] The four panels display the results at z = (2, 3, 4, 5) where 
(F) f = (0.843, 0.601, 0.286, 0.072), corresponding to the values in 
the highest resolution simulation (2.5-400). Note, however, this is 
not necessarily the preferred model due to its small box size. The 
relative difference becomes significantly larger at higher redshifts 
for decreasing mass resolution and box size. 



converges quickly with mass resolution. There is a 0.5 (1.3) 
per cent different between the 20-200 (20-100) and 20-400 
models, giving a good degree convergence with resolution 
for A/ ga s = 1.61 x 10 6 M Q /h. Convergence of (F) with box 
size may be judged by the vertical separation between data 
points with the same gas particle mass. A 20 Mpc/h box pro- 
vides acceptable results at z = 2, with a 0.7 (0.2) percent 
difference between the 20-100 (40-200) and 80-400 models. 
iTVtler et af] l|2009f) . who recently examined the convergence 
of (F) with box size at z — 2 using the Eulerian hydrody- 
namical code Enzo, find similar results. 

However, it is clear that the relative difference becomes 
significantly larger with increasing redshift and decreasing 
mass resolution. By z = 5, there is a 2.1 (8.2) per cent dif- 
ference between {F) in the 5-200 (5-100) and 5-400 models. 
Only a marginal degree of convergence with mass resolution 
is achieved for M gas = 2.01 x 10 5 M /h. The differences due 
to box size also become larger, with a 7.2 (2.1) per cent 
difference between (F) in the 20-100 (40-200) and 80-400 
models. Clearly, the box size and mass resolution require- 
ments for simulating the mean flux of the Lya forest are 
much stricter at higher redshift. 

The second statistic we consider is the ID Lya flux 
power spectrum. This has been extensively used as a probe 
of the primordial matter power spectrum on scales of 0.5 — 
40 Mpc/h at 2 < z < 4, and there are several studies which 
examin e its convergence with resolution an d box size in some 
detail l|McDonaldl 120031 ; IViel et all 120041 '). There has been 
co mparatively little work performed at higher redshifts (but 
see lViel et al.ll2"o"u8h . The upper panels in Fig. [2] display the 
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Figure 2. The upper panels show the ID flux power spectrum 
computed from simulations with fixed box size (lOMpc/h) at z = 
2 and z = 5, while the lower panels display the power spectrum 
from simulations with fixed mass resolution (1.03 X 10 8 M Q /h). 
The spectra have been rescaled to have the same mean flux, indi- 
cated in the upper right of each panel. The difference in the power 
spectrum relative to the models represented by the solid curves 
is shown in the lower third of each panel. T he vertical dotted 
lines bracket the ranee of wavenumbers used bv lViel et all (120041 ) 
to infer the amplitude and shape of the matter power spectrum, 
0.003 < k [skin" 1 ] < 0.03. 



power spectrum of the Lya flux, F = e~ T , at z — 2 and 2 = 5 
computed from the simulations with a box size of 10 Mpc/h. 
The ver tical dotted lines bracket the range of wavenumbers 
used bv lViel et alj ([2004) to infer the amplitude and shape 
of the matter power spectrum at z < 3. Note that we have 
rescaled the synthetic spectra to have the same (F) for this 
comparison. The convergence with mass resolution is again 
significantly poorer at z = 5 and is also scale dependent. 
For the range 0.003 < k [skm -1 ] < 0.03 at z = 2 the 10-200 
(10-100) data is within 1 (3) per cent of the 10-400 model, 
corresponding to a requirement for M gas = 1.61 x 10 Mq/Ii. 
However, at z — 5, the 10-200 (10-100) model deviates from 
10-400 by up to 7 (22) per cent, while the 5-200 (5-100) 
simulations (not shown) deviate from 5-400 by 2 (6) per 
cent. A mass resolution of at least M gas = 2.01 x 10 5 Mq/1i 
is required at z = 5 for marginal convergence. 

The lower panels in Fig. [5] display the effect of box size 
on Pp(fc) for fixed mass resolution. At z = 5, the 20-100 (40- 
200) model is within 10 (8) per cent of the 80-400 model, 
while at z = 2 the 20-100 (40-200) model is within 10 (9) per 
cent of the 80-400 model. A box size of at least 40 Mpc/h 
is therefore required to adequately model the power spec- 
tr um at 2 < z < 5. Sim ilar requirements have been noted 
bv lMcDon ald (2003) and lViel et al. | (120041 1 at z ~ 2 - 3. We 
briefly note we also performed an examination of a third flux 
statistic, the distribution of the Lya flux. We again found 
the simulation requirements to be significantly stricter to- 
wards higher redshifts, with the same box size and resolution 
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Figure 3. The upper panels show the volume weighted gas den- 
sity distribution computed from simulations with fixed box size 
(10 Mpc/h) at z = 2 and z = 5, while the lower panels display the 
density distribution from simulations with fixed mass resolution 
(1.03 X 10 s Mq/1i). The difference in the density distribution rela- 
tive to the models represented by the solid curves is shown in the 
lower third of each panel. The vertical dotted lines in each panel 
bracket the range of optical depth weighted overdensities corre- 
sponding to 95 per cent of all the pixels with 0.05 < F < 0.95 in 
the associated Lya forest spectra. 



requirements as the mean flux. We do not report these re- 
sults in detail here. 

3.2 The gas density distribution 

The explanation for these results is apparent on inspecting 
the gas density distribution in the simulations. The upper 
panels in Fig. display the volume weighted density distri- 
bution per unit log A at z = 2 and z — 5 for models with 
a box size of 10 Mpc/h, while the lower panels display the 
results for fixed mass resolution. The density distributions 
are obtained by interpolating the particle masses, weighted 
by the smoothing kernel, onto a regular grid with a cell size, 
r, indicated in each panel. The vertical dotted lines in each 
panel bracket the range of optical depth weighted overden- 
sities corresponding to 95 per cent of all the pixels with 
0.05 < F < 0.95 in the associated Lya forest spectra. The 
mean flux of the spectra at z — (2, 5) is (F) — (0.865, 0.115). 

At 2 = 2 the 10-200 (10-100) density distribution is 
within 5 (13) per cent of the 10-400 model for —0.5 < 
log A < 2.5, while at z = 5 the 10-200 (10-100) is within 
12 (18) per cent of 10-400, providing at best marginal con- 
vergence. Outwith this density range the distribution has 
not converged at either redshift. In the lower two panels 
the 40-200 (20-100) is within 3 (12) per cent of 80-400 at 
2 = 2, and 7 (18) per cent at z = 5. It is the poor conver- 
gence with mass resolution and box size in the most under- 
dense regions (log A < —0.5) which drives the Lya forest 
results. At 2 = 2, the Lya forest is dominated by gas with 
A > 1. The transmission from underdense regions is always 
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Figure 4. The volume weighted gas density distribution at z = 2 
(left panel) and 2 = 5 (right panel) from the 10-50 (dotted curve) 
and 10-100 (solid curve) models. The dashed curve is computed 
from the 10-100 simulation after doubling the particle smoothing 
lengths, matching those used in the 10-50 model in order to mimic 
lower mass resolution. 



close to the continuum (F = 1) regardless of mass resolu- 
tion, and so under-resolving these regions has little impact 
on the Lya flux statistics. In contrast, underdense regions 
dominate the transmission at z — 5, impacting significantly 
on the convergence of the simulated Lya forest properties. 

This behaviour is a consequence of the spatially adap- 
tive nature of SPH, which provides excellent spatial resolu- 
tion in high density regions but poorer resolution in under- 
dense regions. In Fig. [4] we demonstrate this by comparing 
the density distribution from the 10-50 and 10-100 models 
to a third distribution, again drawn from the 10-100 model. 
The latter is computed by interpolating the particle masses 
onto a grid using the 10-50 model particle smoothing lengths 
(twice the 10-100 values), mimicking particle masses a factor 
of eight larger. The 10-50 distribution at low densities is well 
reproduced by the resmoothed 10-100 distribution, implying 
that differences in the density distribution at log A < —0.5 
are largely a consequence of the intrinsic mass resolution 
limit. However, additional effects such as gas being trans- 
ferred from the low density IGM to previously unres olved 
haloes may also play a small role (|Theuns et al.|[T998h . 

Lastly, we consider fits to the gas density distribution 
which are widely used in analytical models of the Lya forest. 
In Fig. [5] we compare the volume weighted density distribu- 
tion from our 10-400 simulati on to the four paramet e r fit: 
obtained at z = (2, 3, 4) by iMiralda-Escude et af] (|200* 
(hereafter MHR00). The solid curves in Fig. [5] correspond 
to our 10-400 simulation data, while the dashed curves show 
the fits obtained by MHR00. The dot-dashed curve at z = 6 
corresponds t o a m ore recent fi10 to an SPH simulation by 
iPawlik et ail (|2009h (hereafter PSS09), who also use the pa- 
rameterisation suggested by MHR00. 



1 The MHR0 fits are derived from the L10 sim- 
ulation of Miralda-Escud e et al. which uses 
(n m ,tt A ,n b h 2 ,h,a s ,n s ) = (0.4,0.6,0.015,0.65,0.79,0.96), 
with a box size of lOMpc/h box and 288 3 cells. This gives an 
average gas mass per cell of 6 X 10 5 Mq. Note the z = 6 MHR00 
distribution is an extrapolation from the lower redshift fits. 

2 PSS09 use a GADGET-2 simula- 
tion with (fi m , f^Ai f^b^ 2 ) h, ag, n s ) = 
(0.258, 0.742, 0.0228, 0.719, 0.796, 0.963), a box size of 6.25 Mpc/h 
and M gas = 1.8 x 10 s M /h (256 3 gas particles). 
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Figure 5. The volume weighted gas density distribution ex- 
tracted from the 10-400 simulation at z = (2, 3, 4, 6) (solid 
curves). The dotted curves correspond to an eighth order polyno- 
mial fit to the simulation data over the range —1 < log A < 2.5. 
The fits obtained by MHR00 correspond to the dashed curves. 
The dot-dashed curve in the panel at z = 6 also shows the recent 
fit presented by PSS09. The differences in the fits relative to the 
simulation data are shown in the lower third of each panel. 



Although our simulation is in reasonable agreement 
with the MHR00 fits for -0.5 < log A < 1, we confirm 
the claim by PSS09 that the power law tail in the MHR00 
parameterisation, Pv(A) oc A - ' 3 for A 3> 1, provides a 
poor description of the density distribution. This is per- 
haps not too surprising; the MHR00 prescription is based 
on the assumption of a power-law density profile for col- 
lapsed objects and is not obtained directly from the sim- 
ulations. The 10-400 model also differs considerably from 
the MHR00 fits at log A < -0.5, although note that our 
data are not fully converged here. The PSS09 fit is in poor 
agreement with our simulation at z = 6. However, PSS09 
find a similar discrepancy between their fit and simulation 
results, suggesting that the difference between their sim- 
ulated density distribution and our data is actually much 
smaller. Furthermore, we use a very different star formation 
prescription to PSS09 which is designed to optimise Lya for- 
est simulations. We have verified this has little effect on the 
simulated gas density distribution at log A < 2 when com- 
pared to a more sophisticat ed star formation prescription 
l|Springel fe Hernquisti I2003T ). but this choice will produce 
differences in the density distribution at higher densities. 

We conclude that the MHR00 parameterisation is not 
fully adequate for describing Pv from our simulations at 
log A > 1. Consequently, we provide polynomial fits to the 
10-400 density distribution in Table 2 over the range — 1 < 
log A < 2.5 only (dotted curves in Fig. [5}. We deliberately 
avoid parameterising the data, preferring to instead provide 
an accurate representation of the simulations for reference. 
Note, however, that Pv is still only marginally converged 
with resolution for —0.5 < log A < 2.5 and has not fully 
converged with box size. 
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Table 2. Tabulated coefficients for the eighth order polynomial fits to the gas density distribution from our 10-400 simulation, log[Py(x)] = 
^\oi:r l , where x = log A. Py(A) is normalised to unity and APy(A) = P(log A)/ In 10. Note that the fits are made over the range 
— 1 < log A < 2.5 only, and are all within < 5 per cent of the simulation data for —0.5 < log A < 2.5. 
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OO 


«i 


0.1 


a.3 


0,4 


as 


ae 


a 7 


«8 


7.0 


-0.038744 


-1.193136 


-1.209168 


1.480778 


-1.355202 


0.649847 


-0.093190 


-0.024297 


0.006604 


6.5 


-0.045600 


-1.162781 


-1.187344 


1.312205 


-1.225248 


0.647885 


-0.120218 


-0.014682 


0.005503 


6.0 


-0.059247 


-1.148617 


-1.122149 


1.176996 


-1.178106 


0.661365 


-0.105707 


-0.031420 


0.008956 


5.5 


-0.077189 


-1.115055 


-1.049649 


0.913821 


-1.018321 


0.774770 


-0.259515 


0.025978 


0.001650 


5.0 


-0.089459 


-1.083824 


-1.051325 


0.669231 


-0.709611 


0.777568 


-0.414747 


0.101665 


-0.009361 


4.5 


-0.102984 


-1.100786 


-1.064773 


0.677866 


-0.542385 


0.556266 


-0.298218 


0.072713 


-0.006496 


4.0 


-0.132602 


-1.132429 


-0.970648 


0.719472 


-0.552924 


0.418740 


-0.169566 


0.029132 


-0.001301 


3.5 


-0.163770 


-1.134785 


-0.885486 


0.611691 


-0.404073 


0.374816 


-0.230005 


0.070230 


-0.008295 


.3.0 


-0.203625 


-1.166205 


-0.784556 


0.642483 


-0.324878 


0.199731 


-0.147317 


0.059290 


-0.008639 


2.5 


-0.259451 


-1.190049 


-0.614461 


0.661790 


-0.367005 


0.090334 


-0.034647 


0.021446 


-0.004244 


2.0 


-0.325057 


-1.184408 


-0.442799 


0.566567 


-0.347644 


0.074473 


-0.016128 


0.012290 


-0.002695 



4 CONCLUSIONS 

We perform a large set of cosmological SPH simulations to 
explore the effect of mass resolution and box size on two key 
Lycv forest flux sta tistics. As noted in many previous studies 
(see lMeiksinll2007l for a review), we find a mass resolution of 
M gas = 1.61 x 10 6 Mq /h is more than adequate for simulat- 
ing mean Lya flux and power spectrum at z — 2. However, 
towards higher redshift the mass resolution requirement for 
convergence becomes significantly stricter, with a gas par- 
ticle mass at least 8 times smaller required at z = 5. We 
demonstrate this is largely a consequence of the intrinsic 
resolution limit of SPH simulations in low density regions. 
Although a 20Mpc/h box is adequate for simulating the 
mean flux at z — 2, a 40Mpc/h box is required for the 
power spectrum, and this is preferred for both statistics at 
z = 5. We also briefly demonstrate that the MHR00 pa- 
rameterisation for the gas density distribution, although in 
reasonable agreement with our 10-400 model at moderate 
overdensities, provides a poor description of the simulation 
for A > 10. 

Although our results will hold in general for SPH 
Lya forest simulations, in detail they will only be exact for 
the specific models we present. Convergence requirements 
will always depend on the physical process under consider- 
ation, as well as the precision of the observational data to 
which the simulations are compared. We have also assumed 
that the z > 5 forest is dominated by transmission from un- 
derdense regions in the IGM. This picture is correct if the 
UV background is spatially uniform, but if there are large 
fluctuations in the ionising radiation field at z > 5, localised 
patches of highly ionised gas will complicate this interpreta- 
tion somewhat. Lastly, at lower redshifts where the Lya for- 
est transmission is dominated by mildly overdense regions, 
SPH simulations produce results com parable to those f rom 
state-of-the-art adaptive mesh codes |Regan et al.ll2007h . It 
would be very interesting to perform a similar comparison 
at z ~ 5. 
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